home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / hf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  18.0 KB  |  795 lines
  1. /*
  2.  * hf.c
  3.  *
  4.  * Copyright (C) 1989, 1991, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * hf.c,v 4.1 1994/08/09 07:59:14 explorer Exp
  17.  *
  18.  * hf.c,v
  19.  * Revision 4.1  1994/08/09  07:59:14  explorer
  20.  * Bump version to 4.1
  21.  *
  22.  * Revision 1.1.1.1  1994/08/08  04:52:09  explorer
  23.  * Initial import.  This is a prerelease of 4.0.6enh3, or 4.1 possibly.
  24.  *
  25.  * Revision 4.0.1.1  91/09/29  15:44:53  cek
  26.  * patch1: Error messages missing newline.
  27.  * 
  28.  * Revision 4.0  91/07/17  14:38:15  kolb
  29.  * Initial version.
  30.  * 
  31.  */
  32. #ifdef XDRHF
  33. #include    <rpc/rpc.h>
  34. #include    <fcntl.h>
  35. #define        MAXSIZE        ((u_int)1000000)
  36. #endif
  37. #include "geom.h"
  38. #include "hf.h"
  39.  
  40. static Methods *iHfMethods = NULL;
  41. static char hfName[] = "heighfield";
  42.  
  43. static void integrate_grid(), QueueTri();
  44. static int DDA2D(), CheckCell();
  45. static Float intHftri();
  46. static float minalt(), maxalt();
  47.  
  48. typedef struct {
  49.     int stepX, stepY;
  50.     Float tDX, tDY;
  51.     float minz, maxz;
  52.     int outX, outY;
  53.     Vector cp, pDX, pDY;
  54. } Trav2D;
  55.  
  56. unsigned long HFTests, HFHits;
  57.  
  58. #ifdef XDRHF
  59. readbuf(fdp,buf,count)
  60. int *fdp; void *buf; u_int count;
  61. {
  62.     return read(*fdp,buf,count);
  63. }
  64. #endif
  65.  
  66. Hf *
  67. HfCreate(filename)
  68. char *filename;
  69. {
  70.     Hf         *hf;
  71.     float     val, *maxptr, *minptr;
  72.     int     i, j;
  73. #ifdef XDRHF
  74.     static int fp;
  75.     u_int    sqSize;
  76.     XDR        xdrs;
  77.     float    *buf;
  78. #else
  79.     FILE     *fp;
  80. #endif
  81.  
  82. #ifdef XDRHF
  83.     if ((fp = open(filename,O_RDONLY)) < 0) {
  84. #else
  85.     fp = fopen(filename, "r");
  86.     if (fp == (FILE *)NULL) {
  87. #endif
  88.         RLerror(RL_ABORT, "Cannot open heightfield file \"%s\".\n",
  89.                 filename);
  90.         return (Hf *)NULL;
  91.     }
  92.  
  93.     hf = (Hf *)Malloc(sizeof(Hf));
  94.     /*
  95.      * Make the following an option someday.
  96.      */
  97.     hf->BestSize = BESTSIZE;
  98.     /*
  99.      * Store the inverse for faster computation.
  100.      */
  101.     hf->iBestSize = 1. / (float)hf->BestSize;
  102.     /*
  103.      * Get HF size.
  104.      */
  105. #ifdef XDRHF
  106.     xdrrec_create(&xdrs,0,0,&fp,readbuf,NULL);
  107.     xdrs.x_op = XDR_DECODE;
  108.     if (!xdrrec_skiprecord(&xdrs)) {
  109. #else
  110.     if (fread((char *)&hf->size, sizeof(int), 1, fp) == 0) {
  111. #endif
  112.         RLerror(RL_ABORT, "Cannot read height field size.\n");
  113.         return (Hf *)NULL;
  114.     }
  115.         
  116. #ifdef XDRHF
  117.     buf = NULL;
  118.     if (!xdr_array(&xdrs,&buf,&sqSize,MAXSIZE,
  119.                    (u_int)sizeof(float),xdr_float)) {
  120.         RLerror(RL_ABORT, "xdr_array() failed.\n");
  121.         return (Hf *)NULL;
  122.     }
  123.     hf->size = (int)sqrt((double)sqSize);
  124. #endif
  125.     hf->data = (float **)share_malloc(hf->size * sizeof(float *));
  126.     for (i = 0; i < hf->size; i++) {
  127. #ifdef XDRHF
  128.         hf->data[i] = &buf[i*hf->size];
  129. #else
  130.         hf->data[i] = (float *)share_malloc(hf->size * sizeof(float));
  131. #endif
  132.         /*
  133.          * Read in row of HF data.
  134.          */
  135. #ifndef XDRHF
  136.         if (fread((char *)hf->data[i],sizeof(float),hf->size,fp)
  137.             != hf->size) {
  138.             RLerror(RL_ABORT, "Not enough heightfield data.\n");
  139.             return (Hf *)NULL;
  140.         }
  141. #endif
  142.         for (j = 0; j < hf->size; j++) {
  143.             val = hf->data[i][j];
  144.             if (val <= HF_UNSET) {
  145.                 hf->data[i][j] = HF_UNSET;
  146.                 /*
  147.                  * Don't include the point in min/max
  148.                  * calculations.
  149.                  */
  150.                 continue;
  151.             }
  152.             if (val > hf->maxz)
  153.                 hf->maxz = val;
  154.             if (val < hf->minz)
  155.                 hf->minz = val;
  156.         }
  157.     }
  158. #ifdef XDRHF
  159.     (void)close(fp);
  160. #else
  161.     (void)fclose(fp);
  162. #endif
  163.     /*
  164.      * Allocate levels of grid.  hf->levels = log base BestSize of hf->size
  165.      */
  166.     for (i = hf->size, hf->levels = 0; i > hf->BestSize; i /= hf->BestSize,
  167.                 hf->levels++)
  168.             ;
  169.     hf->levels++;
  170.     hf->qsize = CACHESIZE;
  171.     hf->q = (hfTri **)Calloc((unsigned)hf->qsize, sizeof(hfTri *));
  172.     hf->qtail = 0;
  173.  
  174.     hf->lsize = (int *)share_malloc(hf->levels * sizeof(int));
  175.     hf->spacing = (float *)share_malloc(hf->levels * sizeof(float));
  176.     hf->boundsmax = (float ***)share_malloc(hf->levels * sizeof(float **));
  177.     hf->boundsmin = (float ***)share_malloc(hf->levels * sizeof(float **));
  178.  
  179.     hf->spacing[0] = hf->size -1;
  180.     hf->lsize[0] = (int)hf->spacing[0];
  181.     hf->boundsmax[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  182.     hf->boundsmin[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  183.     /*
  184.      * Compute initial bounding boxes
  185.      */
  186.     for (i = 0; i < hf->lsize[0]; i++) {
  187.         hf->boundsmax[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  188.         hf->boundsmin[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  189.         maxptr = hf->boundsmax[0][i];
  190.         minptr = hf->boundsmin[0][i];
  191.         for (j = 0; j < hf->lsize[0]; j++) {
  192.             *maxptr++ = maxalt(i, j, hf->data) + EPSILON;
  193.             *minptr++ = minalt(i, j, hf->data) - EPSILON;
  194.         }
  195.     }
  196.  
  197.     for (i = 1; i < hf->levels; i++) {
  198.         hf->spacing[i] = hf->spacing[i-1] * hf->iBestSize;
  199.         hf->lsize[i] = (int)hf->spacing[i];
  200.         if ((Float)hf->lsize[i] != hf->spacing[i])
  201.             hf->lsize[i]++;
  202.         hf->boundsmax[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  203.         hf->boundsmin[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  204.         for (j = 0; j < hf->lsize[i]; j++) {
  205.             hf->boundsmax[i][j] = (float *)share_malloc(hf->lsize[i] *
  206.                             sizeof(float));
  207.             hf->boundsmin[i][j] = (float *)share_malloc(hf->lsize[i] *
  208.                             sizeof(float));
  209.         }
  210.         integrate_grid(hf, i);
  211.     }
  212.  
  213.     hf->boundbox[LOW][X] = hf->boundbox[LOW][Y] = 0;
  214.     hf->boundbox[HIGH][X] = hf->boundbox[HIGH][Y] = 1;
  215.     hf->boundbox[LOW][Z] = hf->minz;
  216.     hf->boundbox[HIGH][Z] = hf->maxz;
  217.  
  218.     return hf;
  219. }
  220.  
  221. Methods *
  222. HfMethods()
  223. {
  224.     if (iHfMethods == (Methods *)NULL) {
  225.         iHfMethods = MethodsCreate();
  226.         iHfMethods->create = (GeomCreateFunc *)HfCreate;
  227.         iHfMethods->methods = HfMethods;
  228.         iHfMethods->name = HfName;
  229.         iHfMethods->intersect = HfIntersect;
  230.         iHfMethods->normal = HfNormal;
  231.         iHfMethods->uv = HfUV;
  232.         iHfMethods->bounds = HfBounds;
  233.         iHfMethods->stats = HfStats;
  234.         iHfMethods->checkbounds = TRUE;
  235.         iHfMethods->closed = FALSE;
  236.     }
  237.     return iHfMethods;
  238. }
  239.  
  240. /*
  241.  * Intersect ray with height field.
  242.  */
  243. int
  244. HfIntersect(hf, ray, mindist, maxdist)
  245. Hf *hf;
  246. Ray *ray;
  247. Float mindist, *maxdist;
  248. {
  249.     Vector hitpos;
  250.     Float offset;
  251.     Trav2D trav;
  252.  
  253.     HFTests++;
  254.  
  255.     /*
  256.      * Find where we hit the hf cube.
  257.      */
  258.     VecAddScaled(ray->pos, mindist, ray->dir, &hitpos);
  259.     if (OutOfBounds(&hitpos, hf->boundbox)) {
  260.         offset = *maxdist;
  261.         if (!BoundsIntersect(ray, hf->boundbox, mindist, &offset))
  262.             return FALSE;
  263.         hitpos.x = ray->pos.x + ray->dir.x * offset;
  264.         hitpos.y = ray->pos.y + ray->dir.y * offset;
  265.         hitpos.z = ray->pos.z + ray->dir.z * offset;
  266.     }
  267.     /*
  268.      * Find out in which cell "hitpoint" is.
  269.      */
  270.     if (equal(hitpos.x, 1.))
  271.         hitpos.x -= EPSILON;
  272.     if (equal(hitpos.y, 1.))
  273.         hitpos.y -= EPSILON;
  274.  
  275.     if (ray->dir.x < 0.) {
  276.         trav.stepX = trav.outX = -1;
  277.         trav.tDX = -1. / (ray->dir.x * hf->spacing[hf->levels -1]);
  278.     } else if (ray->dir.x > 0.) {
  279.         trav.stepX = 1;
  280.         trav.outX = hf->lsize[hf->levels -1];
  281.         /*
  282.          * (1./size) / ray
  283.          */
  284.         trav.tDX = 1. / (ray->dir.x * hf->spacing[hf->levels -1]);
  285.     }
  286.  
  287.     if (ray->dir.y < 0.) {
  288.         trav.stepY = trav.outY = -1;
  289.         trav.tDY = -1. / (ray->dir.y * hf->spacing[hf->levels -1]);
  290.     } else if (ray->dir.y > 0.) {
  291.         trav.stepY = 1;
  292.         trav.outY = hf->lsize[hf->levels -1];
  293.         trav.tDY = 1. / (ray->dir.y * hf->spacing[hf->levels -1]);
  294.     }
  295.  
  296.     trav.pDX.x = ray->dir.x * trav.tDX;
  297.     trav.pDX.y = ray->dir.y * trav.tDX;
  298.     trav.pDX.z = ray->dir.z * trav.tDX;
  299.     trav.pDY.x = ray->dir.x * trav.tDY;
  300.     trav.pDY.y = ray->dir.y * trav.tDY;
  301.     trav.pDY.z = ray->dir.z * trav.tDY;
  302.  
  303.     trav.cp = hitpos;
  304.     trav.minz = hf->minz;
  305.     trav.maxz = hf->maxz;
  306.     if (DDA2D(hf, &ray->pos, &ray->dir, hf->levels -1, &trav, maxdist)) {
  307.         HFHits++;
  308.         return TRUE;
  309.     }
  310.     return FALSE;
  311. }
  312.  
  313. /*
  314.  * Traverse the grid using a modified DDA algorithm.  If the extent of
  315.  * the ray over a cell intersects the bounding volume defined by the
  316.  * four corners of the cell, either recurse or perform ray/surface
  317.  * intersection test.
  318.  */
  319. static int
  320. DDA2D(hf, pos, ray, level, trav, maxdist)
  321. Hf *hf;
  322. Vector *pos, *ray;
  323. int level;
  324. Trav2D *trav;
  325. Float *maxdist;
  326. {
  327.     int x, y, size, posZ;
  328.     float **boundsmin, **boundsmax, spacing;
  329.     Float tX, tY;
  330.     Trav2D newtrav;
  331.     Vector nxp, nyp;
  332.  
  333.     size = hf->lsize[level];
  334.     spacing = hf->spacing[level];
  335.  
  336.     posZ = (ray->z > 0.);
  337.  
  338.     x = trav->cp.x * hf->spacing[level];
  339.     if (x == size)
  340.         x--;
  341.     y = trav->cp.y * hf->spacing[level];
  342.     if (y == size)
  343.         y--;
  344.     boundsmax = hf->boundsmax[level];
  345.     boundsmin = hf->boundsmin[level];
  346.  
  347.     if (trav->outX > size) trav->outX = size;
  348.     if (trav->outY > size) trav->outY = size;
  349.     if (trav->outX < 0) trav->outX = -1;
  350.     if (trav->outY < 0) trav->outY = -1;
  351.  
  352.     if (ray->x < 0.) {
  353.         tX = (x /spacing - trav->cp.x) / ray->x;
  354.     } else if (ray->x > 0.)
  355.         tX = ((x+1)/spacing - trav->cp.x) / ray->x;
  356.     else
  357.         tX = FAR_AWAY;
  358.     if (ray->y < 0.) {
  359.         tY = (y /spacing - trav->cp.y) / ray->y;
  360.     } else if (ray->y > 0.)
  361.         tY = ((y+1)/spacing - trav->cp.y) / ray->y;
  362.     else
  363.         tY = FAR_AWAY;
  364.  
  365.     nxp.x = trav->cp.x + tX * ray->x;
  366.     nxp.y = trav->cp.y + tX * ray->y;
  367.     nxp.z = trav->cp.z + tX * ray->z;
  368.  
  369.     nyp.x = trav->cp.x + tY * ray->x;
  370.     nyp.y = trav->cp.y + tY * ray->y;
  371.     nyp.z = trav->cp.z + tY * ray->z;
  372.  
  373.     do {
  374.         if (tX < tY) {
  375.             if ((posZ && trav->cp.z <= boundsmax[y][x] &&
  376.                  nxp.z >= boundsmin[y][x]) ||
  377.                 (!posZ && trav->cp.z >= boundsmin[y][x] &&
  378.                  nxp.z <= boundsmax[y][x])) {
  379.                 if (level) {
  380.                     /*
  381.                      * Recurse -- compute constants
  382.                      * needed for next level.
  383.                      * Nicely enough, this just
  384.                      * involves a few multiplications.
  385.                      */
  386.                     newtrav = *trav;
  387.                     newtrav.tDX *= hf->iBestSize;
  388.                     newtrav.tDY *= hf->iBestSize;
  389.                     newtrav.maxz = boundsmax[y][x];
  390.                     newtrav.minz = boundsmin[y][x];
  391.                     if (ray->x < 0.)
  392.                         newtrav.outX=hf->BestSize*x-1;
  393.                     else
  394.                         newtrav.outX=hf->BestSize*(x+1);
  395.                     if (ray->y < 0.)
  396.                         newtrav.outY=hf->BestSize*y-1;
  397.                     else
  398.                         newtrav.outY=hf->BestSize*(y+1);
  399.                     newtrav.pDX.x *= hf->iBestSize;
  400.                     newtrav.pDX.y *= hf->iBestSize;
  401.                     newtrav.pDX.z *= hf->iBestSize;
  402.                     newtrav.pDY.x *= hf->iBestSize;
  403.                     newtrav.pDY.y *= hf->iBestSize;
  404.                     newtrav.pDY.z *= hf->iBestSize;
  405.                     if (DDA2D(hf,pos,ray,level-1,
  406.                         &newtrav, maxdist))
  407.                         return TRUE;
  408.                 } else if (CheckCell(x,y,hf,ray,pos,maxdist))
  409.                     return TRUE;
  410.             }
  411.             x += trav->stepX;        /* Move in X */
  412.             if (*maxdist < tX || x == trav->outX)
  413.                 /* If outside, quit */
  414.                 return FALSE;
  415.             tX += trav->tDX;    /* Update position on ray */
  416.             trav->cp = nxp;        /* cur pos gets next pos */
  417.             nxp.x += trav->pDX.x;    /* Compute next pos */
  418.             nxp.y += trav->pDX.y;
  419.             nxp.z += trav->pDX.z;
  420.         } else {
  421.             if ((posZ && trav->cp.z <= boundsmax[y][x] &&
  422.                  nyp.z >= boundsmin[y][x]) ||
  423.                 (!posZ && trav->cp.z >= boundsmin[y][x] &&
  424.                  nyp.z <= boundsmax[y][x])) {
  425.                 if (level) {
  426.                     /* Recurse */
  427.                     newtrav = *trav;
  428.                     newtrav.tDX *= hf->iBestSize;
  429.                     newtrav.tDY *= hf->iBestSize;
  430.                     newtrav.maxz = boundsmax[y][x];
  431.                     newtrav.minz = boundsmin[y][x];
  432.                     if (ray->x < 0.)
  433.                         newtrav.outX=hf->BestSize*x-1;
  434.                     else
  435.                         newtrav.outX=hf->BestSize*(x+1);
  436.                     if (ray->y < 0.)
  437.                         newtrav.outY=hf->BestSize*y-1;
  438.                     else
  439.                         newtrav.outY=hf->BestSize*(y+1);
  440.                     newtrav.pDX.x *= hf->iBestSize;
  441.                     newtrav.pDX.y *= hf->iBestSize;
  442.                     newtrav.pDX.z *= hf->iBestSize;
  443.                     newtrav.pDY.x *= hf->iBestSize;
  444.                     newtrav.pDY.y *= hf->iBestSize;
  445.                     newtrav.pDY.z *= hf->iBestSize;
  446.                     if (DDA2D(hf,pos,ray,level-1,
  447.                             &newtrav, maxdist))
  448.                         return TRUE;
  449.                 } else if (CheckCell(x,y,hf,ray,pos,maxdist))
  450.                     return TRUE;
  451.             }
  452.             y += trav->stepY;
  453.             if (*maxdist < tY || y == trav->outY)
  454.                 return FALSE;
  455.             tY += trav->tDY;
  456.             trav->cp = nyp;
  457.             nyp.x += trav->pDY.x;
  458.             nyp.y += trav->pDY.y;
  459.             nyp.z += trav->pDY.z;
  460.         }
  461.     } while ((trav->cp.x <= 1. && trav->cp.y <= 1.) &&
  462.          ((posZ && trav->cp.z <= trav->maxz) ||
  463.          (!posZ && trav->cp.z >= trav->minz)));
  464.  
  465.         /*
  466.          * while ((we're inside the horizontal bounding box)
  467.          *        (usually caught by outX & outY, but
  468.          *         it's possible to go "too far" due to
  469.          *         the fact that our levels of grids do
  470.          *         not "nest" exactly if gridsize%BestSize != 0)
  471.          *      and
  472.          *      ((if ray->z is positive and we haven't gone through
  473.          *       the upper bounding plane) or
  474.          *      (if ray->z is negative and we haven't gone through
  475.          *       the lower bounding plane)));
  476.          */
  477.  
  478.     return FALSE;
  479. }
  480.  
  481. /*
  482.  * Check for ray/cell intersection
  483.  */
  484. static int
  485. CheckCell(x, y, hf, ray, pos, maxdist)
  486. int x, y;
  487. Hf *hf;
  488. Vector *ray, *pos;
  489. Float *maxdist;
  490. {
  491.     hfTri *tri1, *tri2;
  492.     Float d1, d2;
  493.     hfTri *CreateHfTriangle();
  494.  
  495.     d1 = d2 = FAR_AWAY;
  496.  
  497.     if (tri1 = CreateHfTriangle(hf, x, y, x+1, y, x, y+1, TRI1))
  498.         d1 = intHftri(ray, pos, tri1);
  499.     if (tri2 = CreateHfTriangle(hf, x+1, y, x+1, y+1, x, y+1, TRI2))
  500.         d2 = intHftri(ray, pos, tri2);
  501.  
  502.     if (d1 == FAR_AWAY && d2 == FAR_AWAY)
  503.         return FALSE;
  504.  
  505.     if (d1 < d2) {
  506.         if (d1 < *maxdist) {
  507.             hf->hittri = *tri1;
  508.             *maxdist = d1;
  509.             return TRUE;
  510.         }
  511.         return FALSE;
  512.     }
  513.  
  514.     if (d2 < *maxdist) {
  515.         hf->hittri = *tri2;
  516.         *maxdist = d2;
  517.         return TRUE;
  518.     }
  519.     return FALSE;
  520. }
  521.  
  522. static hfTri *
  523. CreateHfTriangle(hf, x1, y1, x2, y2, x3, y3, which)
  524. Hf *hf;
  525. int x1, y1, x2, y2, x3, y3, which;
  526. {
  527.     hfTri *tri;
  528.     hfTri *GetQueuedTri();
  529.     Float xid, yid;
  530.     Vector tmp1, tmp2;
  531.  
  532.     /*
  533.      * Don't use triangles with "unset" vertices.
  534.      */
  535.     if (hf->data[y1][x1] == HF_UNSET ||
  536.         hf->data[y2][x2] == HF_UNSET ||
  537.         hf->data[y3][x3] == HF_UNSET)
  538.         return (hfTri *)0;
  539.  
  540.     xid = (Float)x1 / (Float)(hf->size -1);
  541.     yid = (Float)y1 / (Float)(hf->size -1);
  542.  
  543.     if ((tri = GetQueuedTri(hf, xid, yid, which)) != (hfTri *)0)
  544.         return tri;
  545.  
  546.     tri = (hfTri *)Malloc(sizeof(hfTri));
  547.  
  548.     tri->type = which;
  549.     tri->v1.x = xid;
  550.     tri->v1.y = yid;
  551.     tri->v1.z = hf->data[y1][x1];
  552.     tri->v2.x = (Float)x2 / (Float)(hf->size-1);
  553.     tri->v2.y = (Float)y2 / (Float)(hf->size-1);
  554.     tri->v2.z = hf->data[y2][x2];
  555.     tri->v3.x = (Float)x3 / (Float)(hf->size-1);
  556.     tri->v3.y = (Float)y3 / (Float)(hf->size-1);
  557.     tri->v3.z = hf->data[y3][x3];
  558.  
  559.     tmp1.x = tri->v2.x - tri->v1.x;
  560.     tmp1.y = tri->v2.y - tri->v1.y;
  561.     tmp1.z = tri->v2.z - tri->v1.z;
  562.     tmp2.x = tri->v3.x - tri->v1.x;
  563.     tmp2.y = tri->v3.y - tri->v1.y;
  564.     tmp2.z = tri->v3.z - tri->v1.z;
  565.  
  566.     (void)VecNormCross(&tmp1, &tmp2, &tri->norm);
  567.  
  568.     tri->d = -dotp(&tri->v1, &tri->norm);
  569.  
  570.     QueueTri(hf, tri);
  571.     return tri;
  572. }
  573.  
  574. /*
  575.  * Intersect ray with right isoscoles triangle, the hypotenuse of which
  576.  * has slope of -1.
  577.  */
  578. static Float
  579. intHftri(ray, pos, tri)
  580. hfTri *tri;
  581. Vector *pos, *ray;
  582. {
  583.     Float u, v, dist, xpos, ypos;
  584.  
  585.     u = dotp(&tri->norm, pos) + tri->d;
  586.     v = dotp(&tri->norm, ray);
  587.  
  588.     if ((u <= 0. || v > -EPSILON) && (u >= 0. && v < EPSILON))
  589.         return FAR_AWAY;
  590.  
  591.     dist = -u / v;
  592.  
  593.     if (dist < EPSILON)
  594.         return FAR_AWAY;
  595.  
  596.     xpos = pos->x + dist*ray->x;
  597.     ypos = pos->y + dist*ray->y;
  598.  
  599.     if (tri->type == TRI1 && xpos >= tri->v1.x && ypos >= tri->v1.y &&
  600.         xpos + ypos <= tri->v2.x + tri->v2.y)
  601.         return dist;
  602.     if (tri->type == TRI2 && xpos <= tri->v2.x && ypos <= tri->v2.y &&
  603.         xpos + ypos >= tri->v1.x + tri->v1.y)
  604.         return dist;
  605.     return FAR_AWAY;
  606. }
  607.  
  608. /*
  609.  * Compute normal to height field.
  610.  */
  611. /*ARGSUSED*/
  612. int
  613. HfNormal(hf, pos, nrm, gnrm)
  614. Hf *hf;
  615. Vector *pos, *nrm, *gnrm;
  616. {
  617.     *gnrm = *nrm = hf->hittri.norm;
  618.     return FALSE;
  619. }
  620.  
  621. /*ARGSUSED*/
  622. void
  623. HfUV(hf, pos, norm, uv, dpdu, dpdv)
  624. Hf *hf;
  625. Vector *pos, *norm, *dpdu, *dpdv;
  626. Vec2d *uv;
  627. {
  628.     uv->u = pos->x;
  629.     uv->v = pos->y;
  630.     if (dpdu) {
  631.         dpdu->x = 1.;
  632.         dpdu->y = dpdv->z = 0.;
  633.         dpdv->x = dpdv->z = 0.;
  634.         dpdv->y = 1.;
  635.     }
  636. }
  637.  
  638. /*
  639.  * Compute heightfield bounding box.
  640.  */
  641. void
  642. HfBounds(hf, bounds)
  643. Hf *hf;
  644. Float bounds[2][3];
  645. {
  646.     /*
  647.      * By default, height fields are centered at (0.5, 0.5, 0.)
  648.      */
  649.     bounds[LOW][X] = bounds[LOW][Y] = 0;
  650.     bounds[HIGH][X] = bounds[HIGH][Y] = 1;
  651.     bounds[LOW][Z] = hf->minz;
  652.     bounds[HIGH][Z] = hf->maxz;
  653. }
  654.  
  655. /*
  656.  * Build min/max altitude value arrays for the given grid level.
  657.  */
  658. static void
  659. integrate_grid(hf, level)
  660. Hf *hf;
  661. int level;
  662. {
  663.     int i, j, k, l, ii, ji;
  664.     float max_alt, min_alt;
  665.     float **maxinto, **mininto, **frommax, **frommin, *minptr, *maxptr;
  666.     int insize, fromsize, fact;
  667.  
  668.     maxinto = hf->boundsmax[level];
  669.     mininto = hf->boundsmin[level];
  670.     insize = hf->lsize[level];
  671.     frommax = hf->boundsmax[level-1];
  672.     frommin = hf->boundsmin[level-1];
  673.     fact = hf->BestSize;
  674.     fromsize = hf->lsize[level-1];
  675.  
  676.     ii = 0;
  677.  
  678.     for (i = 0; i < insize; i++) {
  679.         ji = 0;
  680.         for (j = 0; j < insize; j++) {
  681.             max_alt = HF_UNSET;
  682.             min_alt = -HF_UNSET;
  683.             for (k = 0; k <= fact; k++) {
  684.                 if (ii+k >= fromsize)
  685.                     continue;
  686.                 maxptr = &frommax[ii+k][ji];
  687.                 minptr = &frommin[ii+k][ji];
  688.                 for (l = 0; l <= fact; l++,maxptr++,minptr++) {
  689.                     if (ji+l >= fromsize)
  690.                         continue;
  691.                     if (*maxptr > max_alt)
  692.                         max_alt = *maxptr;
  693.                     if (*minptr < min_alt)
  694.                         min_alt = *minptr;
  695.                 }
  696.             }
  697.             maxinto[i][j] = max_alt + EPSILON;
  698.             mininto[i][j] = min_alt - EPSILON;
  699.             ji += fact;
  700.         }
  701.         ii += fact;
  702.     }
  703. }
  704.  
  705. /*
  706.  * Place the given triangle in the triangle cache.
  707.  */
  708. static void
  709. QueueTri(hf, tri)
  710. Hf *hf;
  711. hfTri *tri;
  712. {
  713.     if (hf->q[hf->qtail])        /* Free old triangle data */
  714.         free((voidstar)hf->q[hf->qtail]);
  715.     hf->q[hf->qtail] = tri;        /* Put on tail */
  716.     hf->qtail = (hf->qtail + 1) % hf->qsize;    /* Increment tail */
  717. }
  718.  
  719. /*
  720.  * Search list of cached trianges to see if this triangle has been
  721.  * cached.  If so, return a pointer to it.  If not, return null pointer.
  722.  */
  723. static hfTri *
  724. GetQueuedTri(hf, x, y, flag)
  725. Hf *hf;
  726. Float x, y;
  727. int flag;
  728. {
  729.     register int i;
  730.     register hfTri **tmp;
  731.  
  732.     for (i = 0, tmp = hf->q; i < hf->qsize; i++, tmp++) {
  733.         if (*tmp && (*tmp)->v1.x == x && (*tmp)->v1.y == y &&
  734.             (*tmp)->type == flag)
  735.             return *tmp;    /* vertices & flag match, return it */
  736.     }
  737.  
  738.     return (hfTri *)0;
  739. }
  740.  
  741. /*
  742.  * Return maximum height of cell indexed by y,x.  This could be done
  743.  * as a macro, but many C compliers will choke on it.
  744.  */
  745. static float
  746. minalt(y,x,data)
  747. int x, y;
  748. float **data;
  749. {
  750.     float  min_alt;
  751.  
  752.     min_alt = min(data[y][x], data[y+1][x]);
  753.     min_alt = min(min_alt, data[y][x+1]);
  754.     min_alt = min(min_alt, data[y+1][x+1]);
  755.     return min_alt;
  756. }
  757.  
  758. /*
  759.  * Return maximum cell height, as above.
  760.  */
  761. static float
  762. maxalt(y,x,data)
  763. int x, y;
  764. float **data;
  765. {
  766.     float  max_alt;
  767.  
  768.     max_alt = max(data[y][x], data[y+1][x]);
  769.     max_alt = max(max_alt, data[y][x+1]);
  770.     max_alt = max(max_alt, data[y+1][x+1]);
  771.     return max_alt;
  772. }
  773.  
  774. char *
  775. HfName()
  776. {
  777.     return hfName;
  778. }
  779.  
  780. void
  781. HfStats(tests, hits)
  782. unsigned long *tests, *hits;
  783. {
  784.     *tests = HFTests;
  785.     *hits = HFHits;
  786. }
  787.  
  788. void
  789. HfMethodRegister(meth)
  790. UserMethodType meth;
  791. {
  792.     if (iHfMethods)
  793.         iHfMethods->user = meth;
  794. }
  795.